library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(readxl)
library(leaflet)
options(
tigris_class = "sf"
)
Sys.setenv(CENSUS_KEY="10dcd73d7c043e91bac9fb8d3989cbff54b08790")
Here I examine the total number of hospital beds in facilities with the designation “General Acute Care Hospitals” from the California Department of Public Health data, per 1000 members of the population of age 65 and older in each county in California.
First we can look generally at the percent of the population that is of age 65 and older by county in California.
# If not already done, obtain the census data
# acs_vars <-
# listCensusMetadata(
# name = "2018/acs/acs5",
# type = "variables"
# )
# # Save an .rds version on my computer
# setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
# saveRDS(acs_vars, file = "censusData2018_acs_acs5.rds")
# setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")
# obtain the saved census data
setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
acs_vars = readRDS("censusData2018_acs_acs5.rds")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")
# get sex by age for CA counties, modified from the rBasics.Rmd code
CA_county_sexbyage <-
getCensus(
name = "acs/acs5",
vintage = 2018,
region = "county:*",
regionin = "state:06",
vars = "group(B01001)"
) %>%
mutate(
countyval =
paste0(state,county)
) %>%
select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
gather(
key = "variable",
value = "estimate",
-countyval
) %>%
mutate(
label = acs_vars$label[match(variable,acs_vars$name)]
) %>%
select(-variable) %>%
separate(
label,
into = c(NA,NA,"sex","age"),
sep = "!!"
)
# Now find the population that is age 65 and up - also modified from rbasics.Rmd
CA_county_elderly <-
CA_county_sexbyage %>%
filter(!is.na(age)) %>%
mutate(
elderly =
ifelse(
age %in% c(
"65 and 66 years",
"67 to 69 years",
"70 to 74 years",
"75 to 79 years",
"80 to 84 years",
"85 years and over"
),
estimate,
NA
)
) %>%
group_by(countyval) %>%
summarize(
elderly = sum(elderly, na.rm = T),
total_population = sum(estimate, na.rm = T)
) %>%
mutate(
percent_elderly = elderly/total_population
)
# get the counties as sf objects in CA
CA_counties <- counties("CA", cb = F, progress_bar=F)
# map percent elderly
CA_county_elderly %>%
left_join(CA_counties %>% dplyr::select(GEOID), by = c("countyval" = "GEOID")) %>%
st_as_sf() %>%
filter(!is.na(percent_elderly)) %>%
mapview(zcol = "percent_elderly", layer.name = 'Percent ages 65+')
Now let’s combine this data with the hospital data and map the total number of beds in General Acute Care Hospitals per 1000 members of the population ages 65 and above.
# get data on hospital beds
setwd("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Data Library/CHHS")
hospBeds <- read_excel("healthcare_facility_beds.xlsx")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")
# filter by general acute care hospitals
gen_acute_beds <- hospBeds %>% filter(FAC_FDR == 'GENERAL ACUTE CARE HOSPITAL') %>%
group_by(FACNAME, COUNTY_NAME, BED_CAPACITY_TYPE) %>% summarize(gen_acute_beds_hosp_type = sum(BED_CAPACITY, na.rm = T))
# number of beds at each hospital
gen_acute_beds_by_hosp <- gen_acute_beds %>% ungroup() %>% group_by(FACNAME, COUNTY_NAME) %>%
summarize(total_acute_beds_per_hosp = sum(gen_acute_beds_hosp_type, na.rm = T))
# number of beds in each county
gen_acute_beds_county <- gen_acute_beds_by_hosp %>% ungroup() %>%
group_by(COUNTY_NAME) %>%
summarize(totalHosp = n(), total_acute_hosp_beds = sum(total_acute_beds_per_hosp, na.rm = T))
# add county names in upper case to be able to join
CA_counties <- CA_counties %>% mutate(uppercasename = toupper(NAME))
# join county data and beds data
gen_acute_beds_county_w_geom <- CA_counties %>% left_join(gen_acute_beds_county, by = c('uppercasename' = 'COUNTY_NAME')) %>%
mutate_all(~(ifelse(is.na(.), 0, .)))
# join county data and elderly data and get beds per population and per elderly population
gen_acute_beds_county_final <- gen_acute_beds_county_w_geom %>%
left_join(CA_county_elderly, by = c("GEOID" = "countyval")) %>%
mutate(beds_per_1000_elderly = 1000*total_acute_hosp_beds / elderly, beds_per_1000 = 1000*total_acute_hosp_beds / total_population)
# map beds per elderly population
gen_acute_beds_county_final %>%
dplyr::select(beds_per_1000_elderly) %>% st_as_sf() %>% st_set_crs(4269) %>%
mapview(layer.name = "Acute care hospital beds per 1000 residents ages 65+")
Let’s also look at the number of general acute care hospital beds per 1000 members of the total population, just for comparison.
gen_acute_beds_county_final %>%
dplyr::select(beds_per_1000) %>% st_as_sf() %>% st_set_crs(4269) %>%
mapview(layer.name = "Acute care hospital beds per 1000 residents")